<html>
<head>
<title>Blitz++ User's Guide </title>
</head>
<body fgcolor="#27408B" bgcolor="#FFFAF0"  >
<hr>
<ul>
    <li> <a href="blitz05.html">Next chapter</a>
    <li> <a href="blitz03.html">Previous chapter</a>
    <li> <a href="blitz.html">Table of contents</a>
</ul>
<hr>

<a name="l80"></a>
<h1>Chapter 4: Stencils</h1>
<p><a name="arrays-stencils"></a>
    
<!-- BZINDEX stencil objects --><a name="index00382">
<!-- BZINDEX Array!stencils --><a name="index00383">
<p>Blitz++ provides an implementation of stencil objects which
is currently <strong>experimental</strong>.  This means that the exact
details of how they are declared and used may change in future
releases.  Use at your own risk.
<p><br><br><br><table width="100%" border="0" cellpadding=10 align=center><tr><td align="left" bgcolor="#0b6698"><font color="#ffffff" face="Helvetica" size=+5>4.1: Motivation: a nicer notation for stencils</font></td></tr></table><br><a name="l81"></a>

<p>Suppose we wanted to implement the 3-D acoustic wave equation
using finite differencing.  Here is how a single iteration would
look using subarray syntax:
<p><pre>
  Range I(1,N-2), J(1,N-2), K(1,N-2);

  P3(I,J,K) = (2-6*c(I,J,K)) * P2(I,J,K)
    + c(I,J,K)*(P2(I-1,J,K) + P2(I+1,J,K) + P2(I,J-1,K) + P2(I,J+1,K)
    + P2(I,J,K-1) + P2(I,J,K+1)) - P1(I,J,K);
</pre>
<p>This syntax is a bit klunky.  With stencil objects, the
implementation becomes:
<p><pre>
BZ_DECLARE_STENCIL4(acoustic3D_stencil,P1,P2,P3,c)
  P3 = 2 * P2 + c * Laplacian3D(P2) - P1;
BZ_END_STENCIL

  .
  .

  applyStencil(acoustic3D_stencil(), P1, P2, P3, c);
</pre>
<p><br><br><br><table width="100%" border="0" cellpadding=10 align=center><tr><td align="left" bgcolor="#0b6698"><font color="#ffffff" face="Helvetica" size=+5>4.2: Declaring stencil objects</font></td></tr></table><br><a name="l82"></a>

<p><!-- BZINDEX stencil objects!declaring --><a name="index00384">
<p>A stencil declaration may not be inside a function.
It can appear inside a class declaration (in which
case the stencil object is a nested type).
<p>Stencil objects are declared using the macros
<code>BZ_DECLARE_STENCIL1</code>, <code>BZ_DECLARE_STENCIL2</code>, etc.
The number suffix is how many arrays are involved
in the stencil (in the above example, 4 arrays-- P1,P2,P3,c --
are used, so the macro BZ_DECLARE_STENCIL4 is invoked).
<p>The first argument is a name for the stencil object.  
Subsequent arguments are names for the arrays on
which the stencil operates.
<p>After the stencil declaration, the macro <code>BZ_END_STENCIL</code>
must appear (or the macro <code>BZ_END_STENCIL_WITH_SHAPE</code>, described
in the next section).
<p>In between the two macros, you can have multiple assignment
statements, if/else/elseif constructs, function calls,
loops, etc.
<p>Here are some simple examples:
<p><!-- BZINDEX BZ_DECLARE_STENCIL --><a name="index00385">
<p><pre>
BZ_DECLARE_STENCIL2(smooth2D,A,B)
  A = (B(0,0) + B(0,1) + B(0,-1) + B(1,0) + B(-1,0)) / 5.0;
BZ_END_STENCIL

BZ_DECLARE_STENCIL4(acoustic2D,P1,P2,P3,c)
  A = 2 * P2 + c * (-4 * P2(0,0) + P2(0,1) + P2(0,-1) + P2(1,0) + P2(-1,0))
      - P1;
BZ_END_STENCIL

BZ_DECLARE_STENCIL8(prop2D,E1,E2,E3,M1,M2,M3,cE,cM)
  E3 = 2 * E2 + cE * Laplacian2D(E2) - E1;
  M3 = 2 * M2 + cM * Laplacian2D(M2) - M1;
BZ_END_STENCIL

BZ_DECLARE_STENCIL3(smooth2Db,A,B,c)
  if ((c &gt; 0.0) &amp;&amp; (c &lt; 1.0))
    A = c * (B(0,0) + B(0,1) + B(0,-1) + B(1,0) + B(-1,0)) / 5.0
      + (1-c)*B;
  else
    A = 0;
BZ_END_STENCIL
</pre>
<p>Currently, a stencil can take up to 11 array parameters.
<p>You can use the notation <code>A(i,j,k)</code> to read the element at
an offset <code>(i,j,k)</code> from the current element.  If you omit
the parentheses (i.e. as in "A"), then the current element
is read.
<p>You can invoke <em>stencil operators</em> which calculate finite
differences and laplacians.
<p><br><br><br><table width="100%" border="0" cellpadding=10 align=center><tr><td align="left" bgcolor="#0b6698"><font color="#ffffff" face="Helvetica" size=+5>4.3: Automatic determination of stencil extent</font></td></tr></table><br><a name="l83"></a>

<p>In stencil declarations such as
<pre>
BZ_DECLARE_STENCIL2(smooth2D,A,B)
  A = (B(0,0) + B(0,1) + B(0,-1) + B(1,0) + B(-1,0)) / 5.0;
BZ_END_STENCIL
</pre>
Blitz++ will try to automatically determine the spatial extent
of the stencil.  This will usually work for stencils defined
on integer or float arrays.  However, the mechanism does not work
well for complex-valued arrays, or arrays of user-defined types.
If you get a peculiar error when you try to use a stencil, you
probably need to tell Blitz++ the special extent of the stencil
manually.
<p>You do this by ending a stencil declaration with
BZ_END_STENCIL_WITH_SHAPE:
<p><pre>
BZ_DECLARE_STENCIL2(smooth2D,A,B)
  A = (B(0,0) + B(0,1) + B(0,-1) + B(1,0) + B(-1,0)) / 5.0;
BZ_END_STENCIL_WITH_SHAPE(shape(-1,-1),shape(+1,+1))
</pre>
<p>The parameters of this macro are: a TinyVector (constructed
by the shape() function) containing the lower bounds of the
stencil offsets, and a TinyVector containing the upper bounds.
You can determine this by looking at the the terms in the
stencil and finding the minimum and maximum value of each
index:
<p><pre>
        A = (B(0,  0) 
           + B(0, +1)
           + B(0, -1)
           + B(+1, 0)
           + B(-1, 0)) / 5.0;
             --------
  min indices  -1, -1
  max indices  +1, +1
</pre>
<p><br><br><br><table width="100%" border="0" cellpadding=10 align=center><tr><td align="left" bgcolor="#0b6698"><font color="#ffffff" face="Helvetica" size=+5>4.4: Stencil operators</font></td></tr></table><br><a name="l84"></a>

<p><!-- BZINDEX stencil operators --><a name="index00386">
<p>This section lists all the stencil operators provided by Blitz++.
They assume that an array represents evenly spaced data points
separated by a distance of <code>h</code>.  A 2nd-order accurate operator
has error term O(h^2); a 4th-order accurate operator has
error term O(h^4).
<p>All of the stencils have factors associated with them.  For
example, the <code>central12</code> operator is a discrete first derivative
which is 2nd-order accurate.  Its factor is 2h; this means that
to get the first derivative of an array A, you need to use 
<code>central12(A,firstDim)/(2*h)</code>.  Typically when designing
stencils, one factors out all of the <code>h</code> terms for efficiency.
<p>The factor terms always consist of an integer multiplier (often 1)
and a power of <code>h</code>.  For ease of use, all of the operators
listed below are provided in a second ``normalized'' version in which
the integer multiplier is 1.  The normalized versions have
an <code>n</code> appended to the name, for example <code>central12n</code> is
the normalized version of <code>central12</code>, and has factor <code>h</code>
instead of <code>2h</code>.
<p>These operators are defined in <code>blitz/array/stencilops.h</code> if
you wish to see the implementation.
<p><br><br><table width="100%" border="0" cellpadding="3"><tr><td align="left" bgcolor="#1080BF"><font color="#ffffff" face="Helvetica" size=+3>4.4.1: Central differences</font></td></tr></table><a name="l85"></a>

<p><!-- BZINDEX central differences --><a name="index00387">
<p><dl>
<p><p></p><dt><strong>central12(A,dimension)</strong><dd>: 1st derivative, 2nd order accurate.  Factor: 2h
<p><a name="stencils/central12.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-1</td><td>0</td><td>1</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">-1</font></td><td></td><td bgcolor="#000000"><font color="#ffffff">1</font></td></tr>
</table>


<p><p></p><dt><strong>central22(A,dimension)</strong><dd>: 2nd derivative, 2nd order accurate.  Factor: h^2
<p><a name="stencils/central22.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-1</td><td>0</td><td>1</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">1</font></td><td bgcolor="#000060"><font color="#ffffff">-2</font></td><td bgcolor="#000000"><font color="#ffffff">1</font></td></tr>
</table>


<p><p></p><dt><strong>central32(A,dimension)</strong><dd>: 3rd derivative, 2nd order accurate.  Factor: 2h^3
<p><a name="stencils/central32.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-2</td><td>-1</td><td>0</td><td>1</td><td>2</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">-1</font></td><td bgcolor="#000000"><font color="#ffffff">2</font></td><td></td><td bgcolor="#000000"><font color="#ffffff">-2</font></td><td bgcolor="#000000"><font color="#ffffff">1</font></td></tr>
</table>


<p><p></p><dt><strong>central42(A,dimension)</strong><dd>: 4th derivative, 2nd order accurate.  Factor: h^4
<p><a name="stencils/central42.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-2</td><td>-1</td><td>0</td><td>1</td><td>2</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">1</font></td><td bgcolor="#000000"><font color="#ffffff">-4</font></td><td bgcolor="#000060"><font color="#ffffff">6</font></td><td bgcolor="#000000"><font color="#ffffff">-4</font></td><td bgcolor="#000000"><font color="#ffffff">1</font></td></tr>
</table>


<p><p></p><dt><strong>central14(A,dimension)</strong><dd>: 1st derivative, 4th order accurate.  Factor: 12h
<p><a name="stencils/central14.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-2</td><td>-1</td><td>0</td><td>1</td><td>2</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">1</font></td><td bgcolor="#000000"><font color="#ffffff">-8</font></td><td></td><td bgcolor="#000000"><font color="#ffffff">8</font></td><td bgcolor="#000000"><font color="#ffffff">-1</font></td></tr>
</table>


<p><p></p><dt><strong>central24(A,dimension)</strong><dd>: 2nd derivative, 4th order accurate.  Factor: 12h^2
<p><a name="stencils/central24.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-2</td><td>-1</td><td>0</td><td>1</td><td>2</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">-1</font></td><td bgcolor="#000000"><font color="#ffffff">16</font></td><td bgcolor="#000060"><font color="#ffffff">-30</font></td><td bgcolor="#000000"><font color="#ffffff">16</font></td><td bgcolor="#000000"><font color="#ffffff">-1</font></td></tr>
</table>


<p><p></p><dt><strong>central34(A,dimension)</strong><dd>: 3rd derivative, 4th order accurate.  Factor: 8h^3
<p><a name="stencils/central34.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-2</td><td>-1</td><td>0</td><td>1</td><td>2</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">-8</font></td><td bgcolor="#000000"><font color="#ffffff">13</font></td><td></td><td bgcolor="#000000"><font color="#ffffff">-13</font></td><td bgcolor="#000000"><font color="#ffffff">8</font></td></tr>
</table>


<p><p></p><dt><strong>central44(A,dimension)</strong><dd>: 4th derivative, 4th order accurate.  Factor: 6h^4
<p><a name="stencils/central44.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-2</td><td>-1</td><td>0</td><td>1</td><td>2</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">12</font></td><td bgcolor="#000000"><font color="#ffffff">-39</font></td><td bgcolor="#000060"><font color="#ffffff">56</font></td><td bgcolor="#000000"><font color="#ffffff">-39</font></td><td bgcolor="#000000"><font color="#ffffff">12</font></td></tr>
</table>


<p></dl>
<p>Note that the above are available in normalized versions
<code>central12n</code>, <code>central22n</code>, ..., <code>central44n</code> which
have factors of <code>h</code>, <code>h^2</code>, <code>h^3</code>, or <code>h^4</code> as
appropriate.  
<p>These are available in multicomponent versions: for example,
<code>central12(A,component,dimension)</code> gives the central12 operator for the 
specified component (Components are numbered 0, 1, ... N-1).  
<p><br><br><table width="100%" border="0" cellpadding="3"><tr><td align="left" bgcolor="#1080BF"><font color="#ffffff" face="Helvetica" size=+3>4.4.2: Forward differences</font></td></tr></table><a name="l86"></a>

<p><!-- BZINDEX forward differences --><a name="index00388">
<p><dl>
<p><p></p><dt><strong>forward11(A,dimension)</strong><dd>: 1st derivative, 1st order accurate.  Factor: h
<p><a name="stencils/forward11.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>0</td><td>1</td></tr><tr align=right><td></td><td bgcolor="#000060"><font color="#ffffff">-1</font></td><td bgcolor="#000000"><font color="#ffffff">1</font></td></tr>
</table>


<p><p></p><dt><strong>forward21(A,dimension)</strong><dd>: 2nd derivative, 1st order accurate.  Factor: h^2
<p><a name="stencils/forward21.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>0</td><td>1</td><td>2</td></tr><tr align=right><td></td><td bgcolor="#000060"><font color="#ffffff">1</font></td><td bgcolor="#000000"><font color="#ffffff">-2</font></td><td bgcolor="#000000"><font color="#ffffff">1</font></td></tr>
</table>


<p><p></p><dt><strong>forward31(A,dimension)</strong><dd>: 3rd derivative, 1st order accurate.  Factor: h^3
<p><a name="stencils/forward31.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>0</td><td>1</td><td>2</td><td>3</td></tr><tr align=right><td></td><td bgcolor="#000060"><font color="#ffffff">-1</font></td><td bgcolor="#000000"><font color="#ffffff">3</font></td><td bgcolor="#000000"><font color="#ffffff">-3</font></td><td bgcolor="#000000"><font color="#ffffff">1</font></td></tr>
</table>


<p><p></p><dt><strong>forward41(A,dimension)</strong><dd>: 4th derivative, 1st order accurate.  Factor: h^4
<p><a name="stencils/forward41.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>0</td><td>1</td><td>2</td><td>3</td><td>4</td></tr><tr align=right><td></td><td bgcolor="#000060"><font color="#ffffff">1</font></td><td bgcolor="#000000"><font color="#ffffff">-4</font></td><td bgcolor="#000000"><font color="#ffffff">6</font></td><td bgcolor="#000000"><font color="#ffffff">-4</font></td><td bgcolor="#000000"><font color="#ffffff">1</font></td></tr>
</table>


<p><p></p><dt><strong>forward12(A,dimension)</strong><dd>: 1st derivative, 2nd order accurate.  Factor: 2h
<p><a name="stencils/forward12.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>0</td><td>1</td><td>2</td></tr><tr align=right><td></td><td bgcolor="#000060"><font color="#ffffff">-3</font></td><td bgcolor="#000000"><font color="#ffffff">4</font></td><td bgcolor="#000000"><font color="#ffffff">-1</font></td></tr>
</table>


<p><p></p><dt><strong>forward22(A,dimension)</strong><dd>: 2nd derivative, 2nd order accurate.  Factor: h^2
<p><a name="stencils/forward22.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>0</td><td>1</td><td>2</td><td>3</td></tr><tr align=right><td></td><td bgcolor="#000060"><font color="#ffffff">2</font></td><td bgcolor="#000000"><font color="#ffffff">-5</font></td><td bgcolor="#000000"><font color="#ffffff">4</font></td><td bgcolor="#000000"><font color="#ffffff">-1</font></td></tr>
</table>


<p><p></p><dt><strong>forward32(A,dimension)</strong><dd>: 3rd derivative, 2nd order accurate.  Factor: 2h^3
<p><a name="stencils/forward32.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>0</td><td>1</td><td>2</td><td>3</td><td>4</td></tr><tr align=right><td></td><td bgcolor="#000060"><font color="#ffffff">-5</font></td><td bgcolor="#000000"><font color="#ffffff">18</font></td><td bgcolor="#000000"><font color="#ffffff">-24</font></td><td bgcolor="#000000"><font color="#ffffff">14</font></td><td bgcolor="#000000"><font color="#ffffff">-3</font></td></tr>
</table>


<p><p></p><dt><strong>forward42(A,dimension)</strong><dd>: 4th derivative, 2nd order accurate.  Factor: h^4
<p><a name="stencils/forward42.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>0</td><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td></tr><tr align=right><td></td><td bgcolor="#000060"><font color="#ffffff">3</font></td><td bgcolor="#000000"><font color="#ffffff">-14</font></td><td bgcolor="#000000"><font color="#ffffff">26</font></td><td bgcolor="#000000"><font color="#ffffff">-24</font></td><td bgcolor="#000000"><font color="#ffffff">11</font></td><td bgcolor="#000000"><font color="#ffffff">-2</font></td></tr>
</table>


<p></dl>
<p>Note that the above are available in normalized versions
<code>forward11n</code>, <code>forward21n</code>, ..., <code>forward42n</code> which
have factors of <code>h</code>, <code>h^2</code>, <code>h^3</code>, or <code>h^4</code> as
appropriate.  
<p>These are available in multicomponent versions: for example,
<code>forward11(A,component,dimension)</code> gives the forward11 operator for the 
specified component (Components are numbered 0, 1, ... N-1).
<p><br><br><table width="100%" border="0" cellpadding="3"><tr><td align="left" bgcolor="#1080BF"><font color="#ffffff" face="Helvetica" size=+3>4.4.3: Backward differences</font></td></tr></table><a name="l87"></a>

<p><!-- BZINDEX backward differences --><a name="index00389">
<p><dl>
<p><p></p><dt><strong>backward11(A,dimension)</strong><dd>: 1st derivative, 1st order accurate.  Factor: h
<p><a name="stencils/backward11.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-1</td><td>0</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">-1</font></td><td bgcolor="#000060"><font color="#ffffff">1</font></td></tr>
</table>


<p><p></p><dt><strong>backward21(A,dimension)</strong><dd>: 2nd derivative, 1st order accurate.  Factor: h^2
<p><a name="stencils/backward21.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-2</td><td>-1</td><td>0</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">1</font></td><td bgcolor="#000000"><font color="#ffffff">-2</font></td><td bgcolor="#000060"><font color="#ffffff">1</font></td></tr>
</table>


<p><p></p><dt><strong>backward31(A,dimension)</strong><dd>: 3rd derivative, 1st order accurate.  Factor: h^3
<p><a name="stencils/backward31.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-3</td><td>-2</td><td>-1</td><td>0</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">-1</font></td><td bgcolor="#000000"><font color="#ffffff">3</font></td><td bgcolor="#000000"><font color="#ffffff">-3</font></td><td bgcolor="#000060"><font color="#ffffff">1</font></td></tr>
</table>


<p><p></p><dt><strong>backward41(A,dimension)</strong><dd>: 4th derivative, 1st order accurate.  Factor: h^4
<p><a name="stencils/backward41.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-4</td><td>-3</td><td>-2</td><td>-1</td><td>0</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">1</font></td><td bgcolor="#000000"><font color="#ffffff">-4</font></td><td bgcolor="#000000"><font color="#ffffff">6</font></td><td bgcolor="#000000"><font color="#ffffff">-4</font></td><td bgcolor="#000060"><font color="#ffffff">1</font></td></tr>
</table>


<p><p></p><dt><strong>backward12(A,dimension)</strong><dd>: 1st derivative, 2nd order accurate.  Factor: 2h
<p><a name="stencils/backward12.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-2</td><td>-1</td><td>0</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">1</font></td><td bgcolor="#000000"><font color="#ffffff">-4</font></td><td bgcolor="#000060"><font color="#ffffff">3</font></td></tr>
</table>


<p><p></p><dt><strong>backward22(A,dimension)</strong><dd>: 2nd derivative, 2nd order accurate.  Factor: h^2
<p><a name="stencils/backward22.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-3</td><td>-2</td><td>-1</td><td>0</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">-1</font></td><td bgcolor="#000000"><font color="#ffffff">4</font></td><td bgcolor="#000000"><font color="#ffffff">-5</font></td><td bgcolor="#000060"><font color="#ffffff">2</font></td></tr>
</table>


<p><p></p><dt><strong>backward32(A,dimension)</strong><dd>: 3rd derivative, 2nd order accurate.  Factor: 2h^3
<p><a name="stencils/backward32.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-4</td><td>-3</td><td>-2</td><td>-1</td><td>0</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">3</font></td><td bgcolor="#000000"><font color="#ffffff">-14</font></td><td bgcolor="#000000"><font color="#ffffff">24</font></td><td bgcolor="#000000"><font color="#ffffff">-18</font></td><td bgcolor="#000060"><font color="#ffffff">5</font></td></tr>
</table>


<p><p></p><dt><strong>backward42(A,dimension)</strong><dd>: 4th derivative, 2nd order accurate.  Factor: h^4
<p><a name="stencils/backward42.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-5</td><td>-4</td><td>-3</td><td>-2</td><td>-1</td><td>0</td></tr><tr align=right><td></td><td bgcolor="#000000"><font color="#ffffff">-2</font></td><td bgcolor="#000000"><font color="#ffffff">11</font></td><td bgcolor="#000000"><font color="#ffffff">-24</font></td><td bgcolor="#000000"><font color="#ffffff">26</font></td><td bgcolor="#000000"><font color="#ffffff">-14</font></td><td bgcolor="#000060"><font color="#ffffff">3</font></td></tr>
</table>


<p></dl>
<p>Note that the above are available in normalized versions
<code>backward11n</code>, <code>backward21n</code>, ..., <code>backward42n</code> which
have factors of <code>h</code>, <code>h^2</code>, <code>h^3</code>, or <code>h^4</code> as
appropriate.  
<p>These are available in multicomponent versions: for example,
<code>backward42(A,component,dimension)</code> gives the backward42 operator for the 
specified component (Components are numbered 0, 1, ... N-1).
<p><br><br><table width="100%" border="0" cellpadding="3"><tr><td align="left" bgcolor="#1080BF"><font color="#ffffff" face="Helvetica" size=+3>4.4.4: Laplacian operators</font></td></tr></table><a name="l88"></a>

<p><!-- BZINDEX Laplacian operators --><a name="index00390">

<p><dl>
<p><p></p><dt><strong>Laplacian2D(A)</strong><dd>: 2nd order accurate, 2-dimensional laplacian.  Factor: h^2
<p><a name="stencils/Laplacian2D.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-1</td><td>0</td><td>1</td></tr><tr align=right><td>-1</td><td></td><td bgcolor="#000000"><font color="#ffffff">1</font></td><td></td></tr>
<tr align=right><td>0</td><td bgcolor="#000000"><font color="#ffffff">1</font></td><td bgcolor="#000060"><font color="#ffffff">-4</font></td><td bgcolor="#000000"><font color="#ffffff">1</font></td></tr>
<tr align=right><td>1</td><td></td><td bgcolor="#000000"><font color="#ffffff">1</font></td><td></td></tr>
</table>


<p><p></p><dt><strong>Laplacian3D(A)</strong><dd>: 2nd order accurate, 3-dimensional laplacian.  Factor: h^2
<p><p></p><dt><strong>Laplacian2D4(A)</strong><dd>: 4th order accurate, 2-dimensional laplacian.  Factor: 12h^2
<p><a name="stencils/Laplacian2D4.yo"></a>
    
<table cellpadding=2 rules=all>
<tr align=right><td></td><td>-2</td><td>-1</td><td>0</td><td>1</td><td>2</td></tr><tr align=right><td>-2</td><td></td><td></td><td bgcolor="#000000"><font color="#ffffff">-1</font></td><td></td><td></td></tr>
<tr align=right><td>-1</td><td></td><td></td><td bgcolor="#000000"><font color="#ffffff">16</font></td><td></td><td></td></tr>
<tr align=right><td>0</td><td bgcolor="#000000"><font color="#ffffff">-1</font></td><td bgcolor="#000000"><font color="#ffffff">16</font></td><td bgcolor="#000060"><font color="#ffffff">-60</font></td><td bgcolor="#000000"><font color="#ffffff">16</font></td><td bgcolor="#000000"><font color="#ffffff">-1</font></td></tr>
<tr align=right><td>1</td><td></td><td></td><td bgcolor="#000000"><font color="#ffffff">16</font></td><td></td><td></td></tr>
<tr align=right><td>2</td><td></td><td></td><td bgcolor="#000000"><font color="#ffffff">-1</font></td><td></td><td></td></tr>
</table>


<p><p></p><dt><strong>Laplacian3D4(A)</strong><dd>: 4th order accurate, 3-dimensional laplacian.  Factor: 12h^2
</dl>
<p>Note that the above are available in normalized versions
<code>Laplacian2D4n</code>, <code>Laplacian3D4n</code> which have factors
h^2.
<p><br><br><table width="100%" border="0" cellpadding="3"><tr><td align="left" bgcolor="#1080BF"><font color="#ffffff" face="Helvetica" size=+3>4.4.5: Gradient operators</font></td></tr></table><a name="l89"></a>

<p><!-- BZINDEX gradient operators --><a name="index00391">
<p>These return <code>TinyVector</code>s of the appropriate numeric type and length:
<p><dl>
<p><p></p><dt><strong>grad2D(A)</strong><dd>: 2nd order, 2-dimensional gradient (vector of first derivatives),
generated using the central12 operator.  Factor: 2h
<p><p></p><dt><strong>grad2D4(A)</strong><dd>: 4th order, 2-dimensional gradient, using central14 operator.  Factor: 12h
<p><p></p><dt><strong>grad3D(A)</strong><dd>: 2nd order, 3-dimensional gradient, using central12 operator.  Factor: 2h
<p><p></p><dt><strong>grad3D4(A)</strong><dd>: 4th order, 3-dimensional gradient, using central14 operator.  Factor: 12h
<p></dl>
<p>These are available in normalized versions
<code>grad2Dn</code>, <code>grad2D4n</code>, <code>grad3Dn</code> and <code>grad3D4n</code> which
have factors h.
<p><br><br><table width="100%" border="0" cellpadding="3"><tr><td align="left" bgcolor="#1080BF"><font color="#ffffff" face="Helvetica" size=+3>4.4.6: Jacobian operators</font></td></tr></table><a name="l90"></a>

<!-- BZINDEX Jacobian operators --><a name="index00392">
<p>The Jacobian operators are defined over 3D vector fields only
(e.g. <code>Array&lt;TinyVector&lt;double,3&gt;,3&gt;</code>).  They return a
<code>TinyMatrix&lt;T,3,3&gt;</code> where T is the numeric type of the
vector field.
<p><dl>
<p><p></p><dt><strong>Jacobian3D(A)</strong><dd>: 2nd order, 3-dimensional Jacobian
using the central12 operator.  Factor: 2h.
<p><p></p><dt><strong>Jacobian3D4(A)</strong><dd>: 4th order, 3-dimensional Jacobian 
using the central14 operator.  Factor: 12h.
<p></dl>
<p>These are also available in normalized
versions <code>Jacobian3Dn</code> and <code>Jacobain3D4n</code> which
have factors h.
<p><br><br><table width="100%" border="0" cellpadding="3"><tr><td align="left" bgcolor="#1080BF"><font color="#ffffff" face="Helvetica" size=+3>4.4.7: Grad-squared operators</font></td></tr></table><a name="l91"></a>

<p><!-- BZINDEX Grad-squared operators --><a name="index00393">
<p>There are also grad-squared operators, which return <code>TinyVector</code>s of
second derivatives:
<p><dl>
<p><p></p><dt><strong>gradSqr2D(A)</strong><dd>: 2nd order, 2-dimensional grad-squared (vector
of second derivatives), generated using the central22 operator.  Factor: h^2
<p><p></p><dt><strong>gradSqr2D4(A)</strong><dd>: 4th order, 2-dimensional grad-squared, using
central24 operator.  Factor: 12h^2
<p><p></p><dt><strong>gradSqr3D(A)</strong><dd>: 2nd order, 3-dimensional grad-squared,
using the central22 operator.  Factor: h^2
<p><p></p><dt><strong>gradSqr3D4(A)</strong><dd>: 4th order, 3-dimensional grad-squared, using
central24 operator.  Factor: 12h^2
<p></dl>
<p>Note that the above are available in normalized versions
<code>gradSqr2Dn</code>, <code>gradSqr2D4n</code>, <code>gradSqr3Dn</code>, <code>gradSqr3D4n</code>
which have factors h^2.
<p><br><br><table width="100%" border="0" cellpadding="3"><tr><td align="left" bgcolor="#1080BF"><font color="#ffffff" face="Helvetica" size=+3>4.4.8: Curl operators</font></td></tr></table><a name="l92"></a>

<p><!-- BZINDEX curl operator --><a name="index00394">

<p>These curl operators return scalar values:
<p><dl>
<p><p></p><dt><strong>curl(Vx,Vy)</strong><dd>: 2nd order curl operator using the central12 operator.
Factor: 2h
<p><p></p><dt><strong>curl4(Vx,Vy)</strong><dd>: 4th order curl operator using the central14 operator.
Factor: 12h
<p><p></p><dt><strong>curl2D(V)</strong><dd>: 2nd order curl operator on a 2D vector field
(e.g. <code>Array&lt;TinyVector&lt;float,2&gt;,2&gt;</code>), using the central12 operator.
Factor: 2h
<p><p></p><dt><strong>curl2D4(V)</strong><dd>: 4th order curl operator on a 2D vector field,
using the central12 operator.  Factor: 12h
<p></dl>
<p>Available in normalized forms <code>curln</code>, <code>curl4n</code>, <code>curl2Dn</code>,
<code>curl2D4n</code>.
<p>These curl operators return three-dimensional <code>TinyVector</code>s of the 
appropriate numeric type:
<p><dl>
<p><p></p><dt><strong>curl(Vx,Vy,Vz)</strong><dd>: 2nd order curl operator using the central12 operator.  
Factor: 2h
<p><p></p><dt><strong>curl4(Vx,Vy,Vz)</strong><dd>: 4th order curl operator using the central14 operator.  
Factor: 12h
<p><p></p><dt><strong>curl(V)</strong><dd>: 2nd order curl operator on a 3D vector
field (e.g. <code>Array&lt;TinyVector&lt;double,3&gt;,3&gt;</code>, using the
central12 operator.  Factor: 2h
<p><p></p><dt><strong>curl4(V)</strong><dd>: 4th order curl operator on a 3D vector
field, using the central14 operator.  Factor: 12h
<p></dl>
<p>Note that the above are available in normalized versions <code>curln</code> and
<code>curl4n</code>, which have factors of <code>h</code>.
<p><br><br><table width="100%" border="0" cellpadding="3"><tr><td align="left" bgcolor="#1080BF"><font color="#ffffff" face="Helvetica" size=+3>4.4.9: Divergence operators</font></td></tr></table><a name="l93"></a>

<p>The divergence operators return a scalar value.
<!-- BZINDEX divergence operator --><a name="index00395">

<p><dl>
<p><p></p><dt><strong>div(Vx,Vy)</strong><dd>: 2nd order div operator using the central12 operator.
Factor: 2h
<p><p></p><dt><strong>div4(Vx,Vy)</strong><dd>: 4th order div operator using the central14 operator.
Factor: 12h
<p><p></p><dt><strong>div(Vx,Vy,Vz)</strong><dd>: 2nd order div operator using the central12 operator.
Factor: 2h
<p><p></p><dt><strong>div4(Vx,Vy,Vz)</strong><dd>: 4th order div operator using the central14 operator.
Factor: 12h
<p><p></p><dt><strong>div2D(V)</strong><dd>: 2nd order div operator on a 2D vector field,
using the central12 operator.  Factor: 2h
<p><p></p><dt><strong>div2D4(V)</strong><dd>: 2nd order div operator on a 2D vector
field, using the central14 operator.  Factor: 12h
<p><p></p><dt><strong>div3D(V)</strong><dd>: 2nd order div operator on a 3D vector field,
using the central12 operator.  Factor: 2h
<p><p></p><dt><strong>div3D4(V)</strong><dd>: 2nd order div operator on a 3D vector field
using the central14 operator.  Factor: 12h
<p></dl>
<p>These are available in normalized versions
<code>divn</code>, <code>div4n</code>, <code>div2Dn</code>, <code>div2D4n</code>, <code>div3Dn</code>, and
<code>div3D4n</code> which have factors of <code>h</code>.
<p><br><br><table width="100%" border="0" cellpadding="3"><tr><td align="left" bgcolor="#1080BF"><font color="#ffffff" face="Helvetica" size=+3>4.4.10: Mixed partial derivatives</font></td></tr></table><a name="l94"></a>

<p><!-- BZINDEX mixed partial operators --><a name="index00396">
<dl>
<p><p></p><dt><strong>mixed22(A,dim1,dim2)</strong><dd>: 2nd order accurate, 2nd mixed partial derivative.
Factor: 4h^2
<p><p></p><dt><strong>mixed24(A,dim1,dim2)</strong><dd>: 4th order accurate, 2nd mixed partial derivative.
Factor: 144h^2
<p></dl>
<p>There are also normalized versions of the above, <code>mixed22n</code> and
<code>mixed24n</code> which have factors <code>h^2</code>.
<p><br><br><br><table width="100%" border="0" cellpadding=10 align=center><tr><td align="left" bgcolor="#0b6698"><font color="#ffffff" face="Helvetica" size=+5>4.5: Declaring your own stencil operators</font></td></tr></table><br><a name="l95"></a>

<p><!-- BZINDEX stencil operators!declaring your own --><a name="index00397">
<p>You can declare your own stencil operators using the macro
<code>BZ_DECLARE_STENCIL_OPERATOR1</code>.  For example, here is the
declaration of <code>Laplacian2D</code>:
<p><pre>
BZ_DECLARE_STENCIL_OPERATOR1(Laplacian2D, A)
    return -4*A(0,0) + A(-1,0) + A(1,0) + A(0,-1) + A(0,1);
BZ_END_STENCIL_OPERATOR
</pre>
<p>To declare a stencil operator on 3 operands, use 
the macro BZ_DECLARE_STENCIL_OPERATOR3.  Here is
the declaration of <code>div</code>:
<p><pre>
BZ_DECLARE_STENCIL_OPERATOR3(div,vx,vy,vz)
  return central12(vx,firstDim) + central12(vy,secondDim)
    + central12(vz,thirdDim);
BZ_END_STENCIL_OPERATOR
</pre>
<p>The macros aren't magical; they just declare an inline
template function with the names and arguments you specify.
For example, the declaration of <code>div</code> could also
be written
<p><pre>
template&lt;class T&gt;                              
inline typename T::T_numtype div(T&amp; vx, T&amp; vy, T&amp; vz)   
{
  return central12(vx,firstDim) + central12(vy,secondDim)
    + central12(vz,thirdDim);
}
</pre>
<p>The template parameter <code>T</code> is an iterator type for
arrays.
<p>You are encouraged to use the macros when possible, because
it is possible the implementation could be changed in the
future.
<p>To declare a difference operator, use this syntax:
<p><pre>
BZ_DECLARE_DIFF(central12,A) {
  return A.shift(1,dim) - A.shift(-1,dim);
}
</pre>
<p>The method <code>shift(offset,dim)</code> retrieves the element at
<code>offset</code> in dimension <code>dim</code>.
<p>Stencil operator declarations cannot occur inside a function.  If
declared inside a class, they are scoped by the class.
<p><br><br><br><table width="100%" border="0" cellpadding=10 align=center><tr><td align="left" bgcolor="#0b6698"><font color="#ffffff" face="Helvetica" size=+5>4.6: Applying a stencil</font></td></tr></table><br><a name="l96"></a>

<p><!-- BZINDEX stencil objects!applying --><a name="index00398">
<p>The syntax for applying a stencil is:
<p><pre>
  applyStencil(stencilname(),A,B,C...,F);
</pre>
<p>Where <code>stencilname</code> is the name of the stencil, and 
<code>A,B,C,...,F</code> are the arrays on which the stencil
operates.
<p>For examples, see <code>examples/stencil.cpp</code> and <code>examples/stencil2.cpp</code>.
<p>Blitz++ interrogates the stencil object to find out how
large its footprint is.  It only applies the stencil
over the region of the arrays where it won't overrun
the boundaries.
<p>

<p>

<hr>
<ul>
    <li> <a href="blitz05.html">Next chapter</a>
    <li> <a href="blitz03.html">Previous chapter</a>
    <li> <a href="blitz.html">Table of contents</a>
</ul>
<hr>
</body>
</html>
